Monte Carlo study of multiply crosslinked semiflexible polymer networks 
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I. INTRODUCTION 

Networks of semiflexible polymers have become one of the 
focal points in current soft matter research. The reason for 
this interest is twofold: on the one hand, most relevant struc- 
tural biological materials, both intra- and extracellular, share 
the common architecture of crosslinked semiflexible polymer 
networks. Two archetypical examples are the cytoskeleton 
and the extracellular matrix. At the same time, there is a 
wide- spread realization that semiflexible networks represent 
an interesting soft-matter system in their own right, outside of 
any biological context, resulting in a much more fundamen- 
tal interest in the microscopic and geometrical origins of their 
mechanical behavior. 

The mechano-elastic characteristics of networks of semi- 
flexible polymers have been studied to analyze and charac- 
terize different types of these networks, both in vivo [1 ], and 
in vitro Q. The contributions of theory have been many and 
insightful, but analytical progress has typically only been pos- 
sible in certain limiting cases where simplifying assumptions 
may be believed to hold, most notably the assumption of affine 
deformations OH). At the same time, computer simulations 
have been used to study these networks, but they too have had 
to rely on simplifications - either reducing the system to two 
dimensions and limiting to the small- strain regime 1 5 , 6 1 or ig- 
noring the non-linear nature of the constituent filaments 00. 

We believe that the time is right for more realistic numeri- 
cal modeling of these networks that allows for a detailed mi- 
croscopic look at the relations between structure, geometry 
and mechanical properties. To this end, we present a com- 
puter model to simulate these semiflexible polymer networks 
in three dimensions. Networks are considered to consist of fil- 
aments, described as semiflexible polymers. These filaments 
are crosslinked in various locations, which might induce ex- 
tra bending of filaments, thus increasing the free energy of the 
system. We start with a homogeneous, isotropic initial ran- 




FIG. 1 : Schematic presentation of (part of) a semiflexible network, 
in which lines indicate the filaments and dots the crosslinks. The 
section of the filament between crosslinks i and j has length l c jj and 
end-to-end distance rij. 0^/ denotes the angle between two end-to- 
end vectors of neighboring segments along the same filament. 



dom network with a high free energy, and employ a Monte 
Carlo scheme to relax this network. This approach allows 
us to generate realistic three-dimensional networks contain- 
ing hundreds of crosslinks, which are nonetheless well equili- 
brated and thus represent realistic initial conditions for further 
mechanical loading in three dimensions. The methodology to 
generate such networks is the first main result presented in this 
paper, and is described in the first part of the paper (section II). 
In the second part of this paper (section III), we subject 
these networks to shear, and analyze their behavior as a func- 
tion of the network parameters, e.g. the stiffness and the 
length of the filaments. The results of these computer exper- 
iments are compared with experiments to validate our model, 
and yield novel predictions for the mechanical behavior of 
semiflexible networks. 



II. NETWORK GENERATION AND EQUILIBRATION 

We begin our discussion with a detailed look at the gener- 
ation of our semiflexible networks based on single polymer 
energies, and how the Metropolis-Monte Carlo scheme is im- 
plemented and adapted for our specific purposes. 



A. Network Free Energy 

The networks considered in this manuscript consist of fil- 
aments, which are linked by crosslinks / = l...N c . Fig. [T] 
shows a schematic representation of a part of the network, in- 
dicating important parameters of the network and the nota- 
tion used. Each filament is an inextensible semiflexible chain, 
whose energy in the presence of an external force is given by 



Em = 




dt(s) 



ds 




(1) 



where s is the arc length coordinate running along the fila- 
ment, K is the bending stiffness which is related to the per- 
sistence length l p as k = j3~ l l p , with j3 = l/(ki,T), t(s) is 
the (unit) tangent vector along the filament, and / is the ap- 
plied force, directed along the end-to-end vector of the poly- 
mer. The filamentous contribution to the total energy of a net- 
work is the sum of the energies of all filaments. In this work 
we consider inextensible filaments, thus ignoring backbone 
stretching of the filaments, a deformation that is only relevant 
at high forces for most biopolymers. 

A brief note on our nomenclature: our networks consist of 
(multiply) connected filaments. Each of these filaments is par- 
titioned into segments, which begin and end in crosslinks. A 
filament can thus consist of many segments, but is always a 
single mechanical entity, satisfying persistence not just at the 
segment level but also through crosslinks. 

Each crosslink connects segments of two filaments, and our 
networks are therefore strictly tetrafunctional - albeit with the 
possibility of dangling ends which are discarded (we do not 
take steric avoidance into account). Compared to the other 
scales in the network, crosslinks are assumed to be exceed- 
ingly small so that their only action, effectively, is to force a 
binary bond between two distinct filaments, or remote regions 
of the same filaments. 

In our computer simulations, we store a complete list of 
all positions xi of the crosslinks, a complete list of the con- 
tour lengths l Ci tj of the segments between crosslinks i and j, 
and a connectivity table which lists which segments are linked 
by each of the crosslinks. We do not keep track of the spa- 
tial configuration of a segment between two crosslinks. In- 
stead, we use the exact radial distribution function as com- 
puted from Eq. ([T]) [9] to assign to each segment a contour 
length drawn from the radial distribution function computed 
at the segment's end-to-end length and persistence length. In 
this manner, we can already perform an important part of the 
full ensemble sampling in a straightforward manner: different 



assignments of the contour lengths correspond to different re- 
alizations of semiflexible networks with a prescribed spatial 
distribution of crosslinks. The relative likelihood of a given 
distribution of lengths is computed from the free energy of 
the resultant network, which we compute as follows. 

For a given network realization we partition the free energy 
in an internal segment part F2 and an inter-segment part £3 . 
As stated above, the internal degrees of freedom of the seg- 
ments are integrated out. Thus, we express the free energy of 
a segment as a function of the distance between the crosslinks 
(Xij) and the length of the segment (l c ,ij)- If the applied force 
/ in Eq. ([T]) is positive (i.e., stretching the filament), F^ can 
be computed from Eq. ([T]) by employing a semiflexible ana- 
logue of the Marko-Siggia interpolation formula ifTTIl : an ex- 
pression for this is given in the next section. The semiflexible 
WLC force-extension formula is not particularly accurate for 
negative forces, as the filaments quickly assume configura- 
tions with considerable transverse displacements under com- 
pressive loading. The crucial feature of compressive loading, 
however, is that the forces involved are always considerably 
smaller than those encountered for extensional loads - indeed, 
this asymmetry in the force-extension curve is responsible for 
many mechanical features of semiflexible networks. For neg- 
ative forces, we find that the force-extension is adequately de- 
scribed by an exponential approach to the asymptote set by the 
classical Euler buckling force. Integrating the force-extension 
curve yields the following expression for the energy 



9g(r ;7 ) 2 (5+6g(r ;7 )) 
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where gfaj) is the scaled extension given by 
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These equations are not only computationally convenient, 
they also provide an excellent fit to the full, analytical force- 
extension curves as shown in Fig. [2^, where we plot the force 
vs. the scaled extension g{rij). In addition to the single- 
segment force-extension, we also need to keep track of their 
persistence through crosslinks. There is no analytical formula 
for this contribution, and we have therefore simulated many 
individual filaments to obtain a reliable numerical expression 
for this contribution. If the applied force / in Eq. ([T]) is posi- 
tive (i.e., stretching the filament), it turns out that we can cap- 
ture the essential behavior by 



m 



l p°fjk 



(4) 



where L .. and L ., are the contour lengths of the segments 
and 6(jk is the angle between the two end-to-end vectors of the 
segments. Note that this contribution to the total energy is not 
accompanied by an entropic contribution, since it is defined 
by explicit variables in our network. 
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FIG. 2: (Color Online) Validation of our effective Hamiltonian. a) 
Analytic force-extension curve (straight line) (see [10]) vs. interpo- 
lation formula of force-extension (dotted line), where / is the force 
on the segments and g(r) the relative extension with respect to the 
equilibrium configuration, b) The probability distribution P(6) of the 
angle 6 between two segments, for different values of the length l\ , 
resp. I2 of the first and second segment, measured in units of the per- 
sistence length l p . In order of decreasing peak value, the solid lines 
correspond to (l u l 2 ) = (l p /100,l p /50); (7^/100,^/16.7) coincid- 
ing with (Z p /50,Zp/20); (l p /l0,l p /5); and (l p /lOJ p /16.1) coincid- 
ing with (l p /5,lp/2). The dotted lines are our theoretical approxi- 
mations as given in Eq. El. Note that curves for constant l p /(l\ + h) 
fall on top of each other, both in the measured curves and in our 
approximation. 



To assess the quality of the segment-segment energy func- 
tion, we compare the distribution function of this energy with 
simulations. We simulate a single wormlike chain of length 
L w , at a fixed temperature and a persistence length l p , and 
count the probability P w \ c (0) of an angle 6 between the vec- 
tors ~tl w — r^ and f^ — f\ . Here, N is anywhere on the chain. 
Our approximate expression for this probability is P apv (0) ~ 
N(lp/(l Cl -\- Icn)) exp(— J32?3(0)), in which the energy E3 is 
given by Eq. d4b and N(l p /(l Cl -f / C2 )) is a normalization fac- 
tor. This histogram is plotted in Fig.|2]3. Although the corre- 
spondence is not perfect, this formula does reflect the essen- 
tials of the angle distribution, capturing the broadening and 
shift of its peaks. 

In summary, we attribute to a specific network configu- 
ration an energy which is the sum of single-segment ener- 



gies given by Eqs. ([2), plus a sum over all segment-pair 
energies given by Eq. (4]), which runs over all pairs of seg- 
ments belonging to the same filament and meeting in the same 
crosslinks. 



B. Interpolation Formula for the Segment Free Energy 

Eq. ([I]) enables us to derive an analytic approximation for 
the semiflexible force-extension relation. To simplify nota- 
tion, we will pass to dimensionless quantities, rescaling all 
forces by a factor of 1%/k and all lengths by l p /l 2 . Based upon 
Eq. ([I]), we can express the scaled difference between the total 
rescaled length of the polymer (l c ) and the end-to-end length 
at rescaled force <j> (l^) as fll: 



It. 1(1) 
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which gives: 
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At zero force this gives l c — h = 1/6- With this, we can de- 
fine the differential extension at force (f) (i.e, the incremental 
extension compared to that at zero force) as Si = l§ — 1$. We 
use these equations to construct an interpolation formula for 
<j)(rj c jp), which is the direct analogue of the Marko-Siggia 
interpolation for the WLC lfTTIl . Around Iq Eq. ^ gives as a 
first order approximation: 

(j) = 90<5/~. (7) 

In the large force regime we can expand Eq. ([5]) to yield 

(8) 
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Tying the two asymptotes together yields 
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(9) 



which can be integrated once to yield Eq. ^. Fig. [2^ shows 
the comparison between this formula and the exact solution - 
the difference between the two does not exceed 6%. 



C. Network generation 

The task at hand is obviously to determine network configu- 
rations that minimize the free energy thus defined. To this end, 
we use a Monte Carlo minimization scheme: starting from 
an isotropic, random network, we propose random changes in 
topology, each of which is either accepted or rejected accord- 
ing to the Metropolis criterion. 

The initial network is constructed by placing m nodes with 
random coordinates in a cubic periodic cell. To connect these 






(a) 



(b) 



(c) 



FIG. 3: Schematic representation of the three Monte Carlo moves, 
a) Four crosslinks that are connected as shown in the above figure, 
are randomly selected in the network. Bonds AB and CD are bro- 
ken and bonds AC and BD are created, such that the configuration of 
the lower figure is formed after energy relaxation, b) A crosslink is 
randomly chosen at which crosslinks A and B are part of the same 
filament, as are crosslinks C and D (above figure). Now A and C 
become part of the same filament as do B and D. This alters the 
three-crosslink free energy, £3. c) A randomly chosen length (dl) 
is removed from the length of one segment of a filament and trans- 
ferred to a neighboring segment of the same filament, such that the 
configuration of the lower figure is formed after relaxation. 
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FIG. 4: Representation of a generated network. This network con- 
sists of 333 filaments, each on average crosslinked six times. The 
network is periodic in all three directions. Note that the undulations 
of the segments are not represented. 



nodes into a four-fold coordinated network, we proceed iter- 
atively: we begin by identifying three nodes which are close 
to each other and connect these with a loop of three bonds. 
This loop is then extended one bond at a time: we identify a 
node A which is not fully connected and which is closest to an 
existing bond BC, and then replace this existing bond by two 
bonds AB and AC. This process is repeated until all nodes are 
four- fold connected. 

In the resulting fully fourfold-coordinated network, each 
bond is considered to be a segment of a single, long filament. 



This network as a whole can therefore be considered a sin- 
gle, circular filament which is crosslinked to itself at various 
places. We then proceed to minimize the free energy - com- 
puted as detailed before - of this initial network, using the 
standard local minimization method of damped molecular dy- 
namics. 

The initial network will be highly stressed, and in gen- 
eral far removed from a realistic equilibrium configuration. 
Chiefly, this is due to considerable filament bending, with 
intra-filament bends at crosslinks often exceeding 90 degrees. 
As initial large strides towards an optimal configuration will 
proceed along downhill directions related to the release of pre- 
cisely these dominant bending stresses, we first focus on rear- 
ranging the topology of the network, analogous to the contin- 
uous random network approach, pioneered by Wooten, Winer, 
Weaire 1 12] and further extended and optimized as detailed in 
|[T3l . This is realized by a series of Monte Carlo moves that 
alter the topology; these are moves (a) and (b) in Fig. [3] To 
the initial configuration with a topology l! with minimized 
coordinates x f , we assign a free energy F' as obtained from 
Eq. ([4]) plus a quadratic function around the average bond 
distance, to prevent crosslinks from clustering and to tune 
the final network topology. The average bond distance de- 
termines whether the final network will be densely or loosely 
crosslinked. We then change the topology to L" by one of 
the moves, and relax the network with this new topology, re- 
sulting in the new coordinates x" and a new free energy F" . 
Depending on the change in free energy AF = F" — F f , the 
topological change is accepted or rejected, using the Metropo- 
lis algorithm. Note that in this stage, we assume that the free 
energy of a network with minimized crosslinks coordinates 
is representative for the free energy of all networks with the 
same topology, up to some additive constant that is topology- 
independent. 

Once such topology altering moves no longer significantly 
affect the overall energy - this typically happens in configura- 
tions where the bending angle of the filament in each node is 
on average around 20 degrees - contour lengths are attributed 
to the segments. As explained, for a segment AB with end- 
to-end distance r^, the length 1 Ci ab is drawn from the cor- 
responding distribution for the WLC with the desired persis- 
tence length l p . Next, we chop up the single continuous fil- 
ament into many smaller ones, by random deletion of seg- 
ments under the constraint that all crosslinks stay connected, 
up to the point where the desired number of filaments (or, al- 
ternatively, mean filament length) is reached. This network 
is then further equilibrated with the Monte Carlo moves (b) 
and (c) shown in Fig. |3j each of which is now accepted to a 
comparable degree. To avoid computational instabilities for 
floppy filaments we add a short-range repulsive force between 
crosslinks. A typical network generated with this approach is 
shown in Fig. [4] 



III. MECHANICAL RESPONSE OF THE NETWORK 

The ultimate goal is to understand the relationship between 
the structure of a network and its mechanical properties. In the 
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FIG. 5: (color online) a) Master curve of scaled differential modu- 
lus as a function of scaled shear-strain 7/74, where 74 is the strain 
at which the modulus is four times the initial modulus Go- For 
the curves shown, the average number of crosslinks per filament is 
L/l c = 6. The values for the scaled persistence length l p /l c used 
are 15.7, 3.81 and 0.77, of which only the latter is distinguishable at 
large strains (dotted blue line). Besides, we plotted the scaling from 
affine theory, which overlaps with the other curves. The inset shows 
the original strain- stiffness curves, from top to bottom with scaled 
persistence lengths of 15.7, 3.81 and 0.77. Note that we plot all data 
points and draw a curve through them. However, a couple of the data 
points lie well outside the curve (see sect. 2.3). b) Differential mod- 
ulus of the networks as a function of the scaled shear stress a for 18 
networks with varying l p /l c and L/l c . 



following sections, we explore some of the basic mechanical 
properties of our system in an attempt to check whether well- 
known behavior is correctly reproduced, and simultaneously 
to offer a glimpse of the relevant microscopic processes that 
we are now able to study in detail and their role in the overall 
mechanics. 

The behavior of biopolymer networks under strain depends 
on many experimental parameters, such as the concentration 
of biopolymer, the amount of capping proteins and the con- 
centration and characteristics of binding proteins. In our 
simulations, we can reproduce such changes by varying the 
crosslinking length, the persistence length and the average 
number of crosslinks per filament. In this paper we consider 
networks that consist of 10 3 crosslinks connecting 2.10 3 seg- 
ments. Periodic boundary conditions are applied in all three 
directions. We do not take into account the contributions 
of the dangling ends of filaments, nor do we consider the 
excluded volume. Our networks are typically very densely 
crosslinked, which implies that filaments which are close to 



FIG. 6: (color online) a) Differential modulus K as a function of 
lp/lc • The three curves correspond to different strains 7, from bottom 
to top equal to 0.02, 0.10 and 0.20. For all curves shown, L/l c = 6. 
b) K at 7 = as a function of L/l c for various l p /l c (upper line: 
lp/lc = 15.7, middle line: l p /l c = 3.81, bottom line: l p /l c = 0.77). 



each other have a high probability to be crosslinked, and we 
therefore feel it is justified to assume the crosslinking con- 
straints to dominate over the effect of entanglements. Ad- 
ditionally, most naturally occuring biopolymer networks, as 
well as most in vitro biomimetic experiments considered oc- 
cur at fairly low polymeric volume fractions, further reducing 
the importance of excluded volume. The persistence length in 
the networks we use ranges from l p /l c = 1 to l p /l c = 16 and 
the average length of filaments from L/l c = 3 to L/l c = 20. 
The data graphed in this paper is obtained from single, rep- 
resentative networks as large as practically possible to mini- 
mize the effects of finite system size on the observables that 
we measure. 

The experimental techniques used to probe the mechanical 
response are essentially twofold: on the one hand, in vitro net- 
works are often subjected to global shears in commercial rheo- 
metric setups to characterize their macroscopic visco-elastic 
properties (HUHJ. On the other hand, many experiments fo- 
cus rather on the microscopic processes involved by injecting 
small particles (~ \jlm) in the network to monitor the behav- 
ior at the filament scale of the network [Q31Q2). Our compu- 
tational method allows us to work at both levels by direct and 
simultaneous measurement of the overall stiffness as well as 
all individual displacements and forces in the system, to high 
accuracy. 

We model shearing by virtually displacing all crosslink po- 
sitions affinely by small shear-increments of 0.2%. After each 
shear-increment we allow for non-affine relaxation of all indi- 



vidual crosslinks in order to minimize the free energy of the 
network. During this procedure the forces and displacements 
are recorded and can be used for further analysis of the net- 
work response. It's important to note that we allow for full 
relaxation after each strain increment - this would be appro- 
priate for adiabatically slow shears and should therefore be 
compared to the zero-frequency limit in oscillatory rheology, 
as indeed we shall do. 



A. Strain Stiffening 

To characterize our networks, we first consider the differ- 
ential stiffness, K = da/dy during shear. An important and 
characteristic feature of these networks is their highly non- 
linear stiffening behavior under relatively small shear stresses 
lfT4l . As argued in Ref. |4], all experimental curves of the 
modulus of networks of semiflexible polymers collapse for 
small shears on a master curve by scaling the stiffness by the 
initial stiffness (K/Kq) and scaling the shear by its value at 
which the stiffness is four times the initial stiffness (7/74). 
Fig. [5^ shows the scaled strain- stiffness curves of our net- 
works under shear, where we plotted the differential modulus 
K, as a function of shear for different ratio's between l p and l c 
(the average contour length of the segments). We observe the 
same universal scaled stiffening behavior as observed in ex- 
periments. For comparison, we plotted the theoretical curve 
that incorporates the typical force-extension curve of single 
filaments combined with the assumption that the filaments de- 
form affinely, that has shown to represent this same master 
curve. Note that we do not account for rupture and backbone 
stretching of the filaments, which becomes relevant at larger 
shears. The inset shows the original curves, where one clearly 
sees an increase in the initial stiffness as well as a small de- 
crease in the strain at which the networks start to stiffen by 
increasing the stiffness of individual filaments. In our simu- 
lations, l p /l c « 16 is more or less comparable with an actin 
network with an average distance between crosslinks of 1 jlm 
and an average filament length of 6 jlm. Smaller values of 
lp/l c represent networks of filaments with a lower persistence 
length like fibrin or networks that are less dense. 

Another way to compare our results with experiments is 
to look at the scaling in the large strain limit. By superposi- 
tion of a small oscillatory stress on a prestress, the differential 
modulus can be experimentally measured. From these mea- 
surements it is known that K ~ <7 L5 for large stresses a [16]. 
We plotted K/Kq vs. a/cr c , where o c is the critical stress, 
defined as the intersection between the horizontal low-stress 
regime and the high- stress asymptote. As shown in Fig. [5}), all 
our networks show the same characteristic scaling behavior at 
large shears. Combined with the observed stiffening, this indi- 
cates that we capture the essential physics in our model, both 
at small and large shears. 

Eq. ([7]) indicates that the initial stiffness of individual fila- 
ments scales as £o,fil ~ lp/lc- We ex P ect this scaling behav- 
ior to change when the filaments are placed in a network that 
allows non-affine reorientations, as is the case in our simula- 
tions. Since all filament properties scale with l p /l^, we plot 



K vs lp/1%, see Fig-pk The figure shows that Kq ~ (Ip/fy 1 
which emphasizes that for non-affine deformations the persis- 
tence length of the constituent filaments is less important for 
the overall network behavior, as filament reorientations allow 
for an alternative route to comply with the imposed strains. 
For increasing 7, the steepness of the slope increases, which 
strongly correlates with the stiffening of the networks. 

From experiments 031 El and simulations (HJ [HI it is 
known that the average filament length influences the net- 
work response. In cells many capping proteins are active that 
can control the length of the filaments, thus changing the me- 
chanical properties. We measured the initial stiffness Kq as 
a function of the average filament length L and l p . L/l c can 
be considered as the average number of crosslinks on a fila- 
ment, which we can vary while keeping l c constant. As ex- 
pected, Fig. [6J3 shows a decrease in Kq if the average fila- 
ment length decreases. Segments of the same filament influ- 
ence each others displacements, thus restricting the freedom 
to adapt to stresses. Besides, when crosslinks connect two or 
three segments instead of four, these crosslinks are more flex- 
ible to reorient when sheared. Therefore, networks with short 
filaments are softer during shearing. 

Fig. [6}) also shows that the stiffness becomes nearly zero for 
short filaments, a behavior independent of l p . This decrease 
is related to the percolation of the network. When the fila- 
ments become too short, no real network will be formed. In 
that case, shearing will shear the liquid in which the filaments 
are immersed, but the filaments will not be constrained in their 
movement and thus the stiffness will vanish. Please note that 
we employ a specific procedure to remove material from the 
network to generate increasingly sparse networks, which im- 
plies that the filament length at which the modulus vanishes 
cannot be directly related to experiments. The overall trend, 
however, is representative of real networks. 



B. Non-Affine Behavior and Ordering 

While the system deforms in our simulations, we allow for 
non-affine reorientations of the segments. It has been sug- 
gested that such non-affine deformations greatly alter the me- 
chanical response |7], and indeed we find that this is true. 
First, a few words about the definition and measures of non- 
affinity. In general, an applied macroscopic strain maps any 
material point x in the reference space in the network onto a 
new point x' in the target space. The location of the point in 
the target space may be thought of as arising from a combina- 
tion of an affine deformation and a non-affine contribution: 



x' = A(7)x + A(*,7) 



(10) 



where A (7) is the deformation gradient tensor, which for 
the case we consider - three-dimensional simple shear in the 
x\ -direction - is given by 



A( 7 ) 




(11) 



The observation that the non-affine contribution A(x, y) de- 
pends both on the applied strain and the (original) location of 
the point under consideration immediately raises the question 
of what, precisely, it means for a system to be affine. In the 
strictest sense, an affine system may be defined as one obey- 
ing A(jc, y) = for all x, 7. This definition, however, is highly 
restrictive as it does not allow for any non-affine motion at any 
point. For systems that do behave non-affinely to some extent, 
the most general measure of the extent of this non- affinity was 
shown in [ 20 ] to be the non-affinity correlation function 



Aij(x,x f ;Y,Y) 



(Ai(x,y)Aj(xf,y)), 



(12) 



where the average (...) is over all crosslinking points in the 
network. Both its spatial dependence and the strain depen- 
dence are of interest - in a moment we will investigate the 
strain dependent aspects by focusing on the trajectory that sin- 
gle points trace out during a deformation. Following [20], we 
shall measure to this end the equal argument limit of the trace 
of Aij(x,x / ', 7, /) which we shall call simply A(y): 



A( 7 ) = i(|A(x,7) 2 |). 



100 - 



(13) 




Note that this measure need not approach zero at large strains, 
even though one may expect all segments to become aligned 
with the direction of maximal strain in this limit and experi- 
ence, in effect, a purely extensional strain. The reason A (7) 
does not tend to zero lies in the fact that even though the de- 
formation becomes differentially affine, it does not become 
affine in the absolute sense. To focus on this differential affin- 
ity, which we feel is a more appropriate measure of (asymp- 
totic) affinity, we introduce a second measure by considering 
the differential displacement from the initial point X\ to the 
final point Xf before and after a small strain increment dy: 

x f = A(dy)xi + £(*,-, y+dy). (14) 

We use this to define the differential non-affinity measure 



SHY) 



1 



(dy) 2 



|A(*,y) 2 |>. 



(15) 



This measure does go to zero as y becomes very large. Later, 
our simulations will show that we do not expect this limit to 
be attained in experiments as, for realistic parameter values, 
the system will have failed long before. A and 5 A may be 
expressed in terms of each other, and the latter tending to zero 
implies that asymptotically, A should become constant with 
the magnitude of this constant reflecting the overall strength 
of the past non-affinity. 

Ultimately, we are interested to see to what extent non- 
affinity affects the mechanical response. To monitor this influ- 
ence, we perform a shear without relaxation after each strain 
increment, thus obtaining Affine- Fig- [5 snows both ^ a ffi ne , 
which is independent of L/l c , and K for networks with differ- 
ent filament length, all having l p /l c = 1. As can be seen, even 
for long filaments, the difference between affine deformation 
and non-affine deformation is striking, both for he initial mod- 
ulus Kq and for the onset of stiffening. This puts the so-called 



FIG. 7: (color online) Non-affine behavior of the networks during 
shearing, a) The differential modulus K measured during affine 
deformation (upper, dotted line) and during non-affine deformation 
(from bottom to top: L/l c = 5.0, L/l c = 8.0 and L/l c = 20.0). For all 
networks, l p /l c = 0.77. b) Differential non-affinity 8 A as a function 
of strain for different l p /l c . For all curves, L/l c = 6. To relate 8 A 
to other length scales in the system, we plot 8A/r^, where r c is the 
average distance between crosslinks. A value of 1.0 implies that the 
average non-affine displacement is equal to r c if y would be 1.0. The 
inset shows the stiffness vs. the differential non-affinity. 



linear (i.e., small-strain) regime of network elasticity in a new 
perspective: even though the strains are small, there is always 
a finite amount of non-affinity which greatly affects the overall 
small strain response. It is thus crucial to understand the role 
of non-affinity, even at small strains, to predict the network 
modulus. 

Even though the Hamiltonian of the network remains the 
same, the difference between the strain-stiffness curves is 
striking. These differences can only be due to non-affine be- 
havior of the network, as all other determinants - topology, 
filament length, density and persistence length, are identical. 
There has been some debate whether the origin of stiffening is 
ultimately entropic or mechanical, but our results suggest that 
rather, we should focus our attention on the degree of non- 
affinity which acts to delay and attenuate the stiffening. 

To see whether our systems tend to affinity at the largest 
strains, we measure the differential non-affinity 5 A as a func- 
tion of the applied macroscopic strain. To relate 8 A to other 
length scales in the system, we plot <5A/r;r, where r c is the av- 
erage distance between crosslinks. Fig. [7]) shows a strong in- 
crease in 5 A with increasing strain. From the inset of Fig. [T]}, 
a strong correlation between the stiffening of the networks and 
the amount of non-affine displacements is revealed. Appar- 
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FIG. 8: (color online) Non-affine behavior of the networks during 
shearing, a) The scaled overall non- affinity A /r* at y = as a func- 
tion of L/l c for persistence lengths l p /l c = 0.77 (•) and l p /l c = 15.7 
(*). Curves are drawn as a guide to the eye. b) Ordering during 
shearing, after subtracting the value of CO at y = 0. The solid line 
indicates the ordering due to affine shearing. The dotted lines indi- 
cate the ordering during shearing in networks with (from bottom to 
top) L/l c = 4.0, L/l c = 6.0 and L/l c = 12.0. c) Distribution of an- 
gles with respect to the x-axis, both for an affine deformation (dotted 
bars) and a non-affine deformation (solid bars) at y = 0.7. The curve 
is the analytic expression for the angular distribution of an initially 
isotropic material at shear y= 0.7. The inset shows the sheared box, 
in which two connected crosslinks are indicated by dots; their end- 
to-end vector makes an angle with the i-axis. 



ently, to prevent the extreme extension that filaments would 
experience at high strains in an affine setting, the network 
shows a strong non-affine reorientation. As indicated before, 
we expect that for high shear all filaments will be aligned in 
the direction of shear and deform purely by stretching. Since 
stretching is an affine deformation, we expect 5 A to ultimately 
tend to zero at large strain. This figure and the inset make it 
clear that this asymptotically (differentially) affine regime is 



FIG. 9: (color online) Correlation 8N(r/r c ) of non-affine behavior, 
as a function of distance r/r c ; a) for strains y = 0.0 (K = K$\ bottom 
curve), y = 0.2 (K = 4.0Kq), y = 0.25 (K = 1 .6Kq), y = 0.29 (K = 
25K ) and y = 0.33 (K = 260K ; top curve), b) for networks with 
different l p /l c (solid line: l p /l c = 0.77, dash-dotted line: l p /l c = 
3.81 and dashed line: l p /l c = 15.7); for all networks, L/l c = 6. c) 
for networks with different L/l c , ranging from 4 to 20. All curves 
collapse, except the dash-dotted (brown) curve with L/l c = 4. For all 
networks, l p /l c = 15.7. 



never actually attained, and non-affinity will continue to fea- 
ture prominently all the way up to the point of failure. 

Shorter filaments, or filaments that are less densely 
crosslinked, are less constrained in their motions which 
should, in principle, allow for greater non-affine motions. To 
verify whether indeed this is the case, we evaluate the non- 
affinity as a function of filament length. To better compare to 
existing experiments, we shall use A, the overall non-affinity 
parameter, instead of the differential measure 5 A. Fig. [8^ in- 
deed shows a pronounced increase in non-affinity as the length 
decreases. This is in agreement with experiments on f-actin, 
which also show an increase in A(0) for decreasing filament 
length ifTTl . Translating A(0) to real distances gives for actin 
with an average l c = \jlm values between 2 — 6/lm 2 , which 



is close to the values between 2 — 10/lm 2 reported for exper- 
iments 1 17]. Besides the length dependence, we also observe 
a dependency of l p on the non-affinity: networks of stiff er fil- 
aments behave more non-affinely. This suggests that asymp- 
totically, we recover the classical picture of rubber elasticity 
(which does very well for long, flexible polymers but rather 
poorly for semiflexible systems): as the persistence length de- 
creases, the polymer configurations become increasingly ran- 
dom (i.e., Gaussian) which is accompanied by a decrease in 
the non-affinity. This is precisely the rubber limit: Gaussian 
polymers deforming affinely. Note that these non-affine de- 
formations are the sole possible origin of the large difference 
in stiffness between affine and non-affine deformations shown 
above. Thus, even though the magnitude of the non-affine de- 
formations is small, they do have an important effect on the 
network response. This is a striking example of the value of 
simulations in this field: microscopic structure and motion are 
of crucial importance to properly understand the macroscopic 
behavior. 

Apparently, the macroscopic result of the microscopic non- 
affinity is to lower the overall stiffness of the system. This 
suggests an interesting question: if the filaments do not go 
to their affine positions, where do they go? To begin to 
answer this, we consider the orientational order of our net- 
works and compute the nematic order parameter &>, defined 
as co = (3 cos 2 — l)/2. Here the average is taken over all 
vectors connecting crosslinks that are connected by segments 
of filaments and is the angle between such a vector and the 
average orientation. An isotropic network has CO = 0, while a 
fully ordered network has co = 1 . Even when we shear a net- 
work affinely, the order will increase from zero to one. To ap- 
preciate the effect non-affinity has, we should therefore com- 
pare to the affine ordering. This affine ordering is represented 
by the solid line in Fig. [8J3. The dotted lines show the effect 
of non-affine reorientations on the ordering of the network, 
for different filament lengths. Interestingly, non-affine reori- 
entations tend to increase the order in the network, a behavior 
independent of l p . 

To get insight in the direction of the ordering, we plot the 
distribution of the angle (j) of the end-to-end vectors of seg- 
ments with respect to the x-axis at y = 0.7, as shown in Fig. [SJ;. 
By comparing the distribution in a non-affine network defor- 
mation (straight bars) with the distribution of an affine net- 
work deformation (dotted bars) we see that the non-affinity 
increases the number of segments oriented at a small angle. To 
appreciate the differences in the two distributions, we plot the 
analytic expression for the distribution of an initial isotropic 
medium that is sheared affinely. Interestingly, the maximum 
of P((j>) coincides with the maximal extensional strain expe- 
rienced as a function of angle. As the figure clearly shows, 
the additional ordering is in the direction of maximal exten- 
sional strain. This might seem counterintuitive - the order ap- 
pears to be increasing in the direction of increasing filament 
extensional strain, which would be highly unfavorable from 
an energetic point of view. However, one should keep in mind 
that non-affine motions are not purely rotational: they may 
encompass additional and simultaneous overall shifts and ex- 
tensional/compressional components. It would be most inter- 



esting to see if this increased order is also observed in experi- 
ments. Our simulations suggest that systems containing long 
filaments are the best place to look for this effect, even though 
these tend to display lower overall non-affinity. 

So far, we have considered only the non-affine motion of 
single points. The non-affinity correlation function is not only 
a function of strain, it may also be evaluated for spatially sep- 
arated points x and x' . To this end, we consider SN(r) = 
((r— r a ff) 2 ) r /<5y 2 , where r is the actual vector between two 
crosslinks and r a ff is the vector between the crosslinks if they 
would have moved affine during 8y. The average, now, runs 
over all pairs of crosslinks whose separation is r. 

As explained in ref. ifTTlL there are two limiting cases in 
the behavior of 8N(r). If filaments would be stiff rods, the 
only way to adapt to strain would be by rotating the whole fil- 
ament. In that case, doubling r would double r — r a ff and thus 
8N(r) ~ r 2 . In the other limiting case, segments along a fila- 
ment behave totally uncorrelated, leading to 8N(r) ~ r°. The 
latter is also the limit for r — > oo. However, the net effect of 
correlated motion of segments along a filament will be highly 
sensitive to the actual network configuration. 

Fig. [9^ shows 8N(r) of a network at different strains. Note 
that the larger scatter for small values of r is due to the smaller 
number of pairs of crosslinks. As can be seen, 8N(r) varies 
with strain: low strains give a low initial value of 8N(r) 
and a steep increase while high strains show just the oppo- 
site. This behavior is strongly correlated to the stiffening be- 
havior shown in Fig. [5^ (upper line). The observed strain- 
dependence of 8N(r) is indiscernible when normalizing with 
respect to y rather than of 5y, which might explain why ex- 
periments report no strain-dependence ifTTl . 

We observe a small but systematic dependence on l p , as 
shown in Fig. [9J3. Interestingly, thus far we hardly observe 
any length dependence of the correlation in non-affine behav- 
ior, which is shown in Fig.[9fc for networks with l p /l c = 15.7. 
For filaments with l p > l c , one would naively expect that the 
behavior of segments along a filament will be much more cor- 
related than the behavior of segments belonging to different 
filaments. Thus, one might expect to find increasing spatial 
correlations for systems composed of larger filaments. That 
we do not see this behavior suggests that it is approximately 
balanced by another effect: larger filaments have more links 
to the rest of the system and are therefore more constrained. 
While the individual segments along a single filament would 
like to line up, they become increasingly unable to do so. In- 
terestingly, the first experiments to measure 8N(r) do show a 
length dependence [17). We cannot rule out that we will see 
this behavior at larger system sizes, but for now are unable to 
reproduce it. 



C. Collective Rearrangements 

A closer look at Fig. [5^ reveals some outliers in the K vs. 
y curve. These discontinuities in K are accompanied by an 
increase in A. This is not a glitch, but rather reflects an in- 
teresting microscopic aspect of our networks. To understand 
this behavior we look at the displacements of individual seg- 
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FIG. 10: Illustration of a collective reorientation in one of our net- 
works with l p /l c = 15.7 and L/l c = 6 during deformation at respec- 
tively 7 = 0.268, 7 = 0.270 and y = 0.272. The thickness of the 
segment indicates the size of its displacement. The modulus K be- 
longing to the deformation of this network is shown as the upper 
curve in the inset in Fig. [5^. 



ments during shearing. Fig. 10 shows a network in which the 
thickness of the segments indicates their displacement during 
a strain increment of 0.2%. Here we see what happens: during 
such a strain increment, a significant fraction of the segments 
has a relatively large incremental displacement in comparison 
with the average displacement of segments during shear incre- 
ments. 

The noteworthy feature is not so much that there are large 
displacements, but rather that these displacements are local- 
ized and occur in correlated fashion. This is reminiscent of 
the behavior of so-called collectively rearranging regions, ob- 
served in simulation and experiment in glassy systems and 
colloidal suspensions. These events are rare over the time 
courses that we have simulated, but may turn out to play an 
important role in the long-time behavior of these materials. 
It would be most interesting to check whether these events 
are also seen in experiments - while these may not be able to 



resolve the blip in K they might be able to register the accom- 
panying peaks in A. The weight in determining A of a reori- 
entation of a certain size decreases with increasing 7, since A 
measures the total non-affinity relative to the total shear. This 
implies that for small shears, reorientations might induce huge 
peaks in A, while these peaks are absent for larger shears even 
though the reorientations are still present. 



IV. CONCLUSIONS 

We have presented a new method to generate and deform 
3D networks of biopolymer filaments. By an adequate choice 
of energies both the entropic stiffness of individual segments 
as well as the persistence of filaments through crosslinks can 
be taken into account. By a Monte Carlo thermalization the 
networks find a local minimum, without further interference 
from our side. 

This method enables us to relate the macroscopic network 
response to microscopic behavior of individual segments and 
crosslinks, both at small and large strains. Although a quan- 
titative comparison between experiments and our simulations 
is hard to obtain, the first results from these simulations agree 
well with experiments. Both the stiffening during shearing 
and the length-dependency of the non-affinity are as expected 
and fit well into the general framework of the behavior of 
semiflexible polymers. Besides, the stress-dependence of the 
stiffness for large shears is the same as experiments have 
shown. This confirms that our model captures the right fea- 
tures that decide the network behavior. 

Our model proves an excellent tool to compare affine defor- 
mations with deformations that allow for non-affine displace- 
ments. We have shown that non-affine displacements have a 
large influence on the stiffness of a network and the onset of 
stiffening. This accounts for the important role of filament 
length. Besides, the accuracy of analysis of the behavior of 
the filaments during deformation reveals some surprising re- 
sults that are hard to obtain by experimental analysis. Thus 
far unobserved, the non-affinity increases the order in the net- 
works. 

Thus far we have only considered networks of a single type 
of filaments. Both in cells and in the extracellular matrix, 
the important load-bearing biopolymer networks are made up 
of different kinds of filaments: vessel walls are composed of 
collagen and elastin, and the cytoskeleton too is a composite 
system containing f-actin, intermediate filaments and micro- 
tubules. This method is a promising tool to explore the be- 
havior of such composite networks under strain, and we are 
currently exploring their properties. 
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